function varargout = XglomPipeline(varargin)
% XGLOMPIPELINE MATLAB code for XglomPipeline.fig
%      XGLOMPIPELINE, by itself, creates a new XGLOMPIPELINE or raises the existing
%      singleton*.
%
%      H = XGLOMPIPELINE returns the handle to a new XGLOMPIPELINE or the handle to
%      the existing singleton*.
%
%      XGLOMPIPELINE('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in XGLOMPIPELINE.M with the given input arguments.
%
%      XGLOMPIPELINE('Property','Value',...) creates a new XGLOMPIPELINE or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before XglomPipeline_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to XglomPipeline_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help XglomPipeline

% Last Modified by GUIDE v2.5 07-Feb-2013 15:04:23


%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
%
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.


% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @XglomPipeline_OpeningFcn, ...
    'gui_OutputFcn',  @XglomPipeline_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    
    if strfind(varargin{1},'1008')==1
        filename=varargin{1};
    else
        gui_State.gui_Callback = str2func(varargin{1});
        disp(varargin{1});
    end
else
    varargin{1} = '1008.2.18.3.1.1.1';
end


if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT



% --- Executes just before XglomPipeline is made visible.
function XglomPipeline_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to XglomPipeline (see VARARGIN)

handles.Xglom_version= getenv('XGLOM_VERSION');
if isempty(handles.Xglom_version)
    [status handles.Xglom_version]=system('awk -F''='' ''/XGLOM_VERSION=/ {print $2}'' xglom.sh| tr -d ''\n''')
end

addpath(fullfile(pwd,'plugins','progressbar'));

% Choose default command line output for XglomPipeline
handles.output = hObject;
handles.Selection = 1;
handles.SelectAll = 0;
handles.Datatype='MR'; % default 'MR', others 'CT', or 'VM'
handles.Img=[];
handles.debug=0;
guidata(hObject, handles);


optargin = size(varargin,2);
stdargin = nargin - optargin;


n = 1;
while n <= numel(varargin)
    switch varargin{n}
        case 'IMG'
            if n == numel(varargin)
                error 'Clim vector not supplied'
            end
            if ~isnumeric(varargin{n+1}) || ndims(varargin{n+1}) ~= 3
                error 'Invalid image volume '
            end
            handles.Img = single(varargin{n+1});
            n = n + 2;
        case 'DATATYPE'
            if n == numel(varargin)
                error 'DATATYPE str not supplied'
            end
            if ~ischar(varargin{n+1})
                error 'Invalid DATATYPE string'
            end
            handles.Datatype =varargin{n+1};
            n = n + 2;
        case 'DARISID'
            if n == numel(varargin)
                error 'DaRIS id not supplied'
            end
            if ~ischar(varargin{n+1})
                ShowStatus('XglomRunPipeline:  must have full Daris ID string.',handles)
                error 'Invalid DaRIS id string'
            end
            
            handles.DatasetID=varargin{n+1};
            n = n + 2;
    end
end

ShowStatus(['XglomRunPipeline: Ready to run pipline on dataset DaRIS ID: ' handles.DatasetID], handles);

table = get(handles.uitable1,'Data');

table{1,1}= char(handles.DatasetID);

set(handles.uitable1,'Data',table);

if ~isempty(handles.Img)
    midaxis = floor(size(handles.Img,1)/2);
    iH = imagesc(single(squeeze(handles.Img(midaxis,:,:))), 'parent', handles.pipelineaxes);
    set(iH, 'hittest', 'off');colormap(gray);
    axis(handles.pipelineaxes, 'equal','tight');
    set(handles.pipelineaxes, ...
        'box'             , 'on', ...
        'xtickmode'       , 'auto', ...
        'ytickmode'           , 'auto', ...
        'interruptible'   , 'off', ...
        'cameraviewanglemode' , 'auto', ...
        'dataaspectratiomode' , 'auto', ...
        'busyaction'      , 'queue');
end

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes XglomPipeline wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = XglomPipeline_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;



% --- Executes during object creation, after setting all properties.
function listbox1_CreateFcn(hObject, eventdata, handles) %#ok<*INUSD,*DEFNU>
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- Executes on button press in pushbutton1. PIPELINE
function pushbutton1_Callback(hObject, eventdata, handles)
% PUSHBUTTON1_CALLBACK calls the pipeline processing function
% reguired for the data set.
% TODO: Fix bug for different data types, e.g. CT

ButtonString = get(hObject,'String')
switch ButtonString
    case 'Next Step'
        uiresume(gcbf)
        
    case 'Begin'
        set(hObject,'String','Next Step')
        switch handles.Datatype
            case 'MR'
                %   info=GetDICOMInfo(handles.DatasetID)
                
                %if isempty(strfind(info.ImageType,'ORIGINAL'))
                % Pre-filtered MR image of original Agilent 9.4 MR image
                %     XglomRunPipeline_version1(hObject, eventdata, handles)
                % else
                % Primary image generated by Agilent 9.4 MR console
                XglomRunPipeline_version2(hObject, eventdata, handles)
                % end
            case 'CT'
                XglomRunPipeline_CTversion(hObject, eventdata, handles)
        end
        
end


function XglomRunPipeline_version1(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
ShowStatus('XglomRunPipeline: Running Pipeline ... ',handles);
% Selection = handles.Selection;
% DatasetID='';
% if isempty(Selection) %|| length(Selection) > 3
%     disp(' No dataset selected.')
%     return
% else
%     table = get(handles.uitable1,'Data');
%     DatasetID = table{Selection,1};
%     if ~isempty(DatasetID)
%         img = GetDICOMImage(DatasetID);
%     else
%         display ' No dataset'
%         return
%     end
% end
%
% ax = get(handles.axes1)
% imagesc(squeeze(img(256,:,:)),'Parent',ax);colormap(gray); axis image; axis on;

tstart=tic;
img = GetDICOMImage(handles.DatasetID);
%imagesc3(img)
gui_active(1);      % will add an abort button
hProgress=figure;
ShowStatus('XglomRunPipeline: Phase 1: Image Processing Procedure (approx. 4 mins)',handles)
hProgress           = progressbar([],0,' Phase 1: Image Processing Procedure (approx. 4 mins)' );
totalprocessing = 1140; %estimate 240 secs for imfill, 18 secs each for bwlabel

%% Phase 1: Image filling
telapsed = toc(tstart);
imf=single(imfill(img,'holes')); %default:conn=6
depth=imf-single(img);
%load ./Cache/1008.2.18.2.1.1.1/depth_mask_01.mat
ShowStatus('XglomRunPipeline: Phase 1: Normalising Processed Image',handles);
hProgress = progressbar( hProgress, 230/totalprocessing, 'Phase 1: Normalising Processed Image');
[mask newimg] =KidneyMask(img,1,497,503);
save(['./Cache/' handles.DatasetID '/depth_mask_01.mat'],'depth','mask');


MatrixSize=[512 512]; %Size of the Depth Dataset

depth = depth + abs(min(depth(:)));
depth = depth./max(depth(:));
%imagesc3(depth)

%% Phase 2: Thresholding
ShowStatus('XglomRunPipeline: Phase 2: Performing Threshold Statistics (x 50)')
hProgress = progressbar( hProgress, 240/totalprocessing, 'Phase 2: Performing Threshold Statistics (x 50)');
processingcount=240;


telapsed = toc(tstart);

ShowStatus('XglomRunPipeline: Running Statistical tests ...', handles);
% addpath('/gpfs/M2Home/projects/Monash016/RatKidney/Agilent')

DepthDataset = depth;
ThresholdRange= 0.01:0.01:0.5;
stored_hist=zeros(50,1000);
stored_diams = zeros(50,141);
counter=1;
for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %%
    %Thresholding
    Holes=DepthDataset>ThresholdRange(counter);
    %Discard single voxels
    ConArray=conndef(3,'minimal');
    ConArray(2,2,2)=0;
    Holes=imdilate(Holes,ConArray).*Holes;
    Holes=single(Holes);
    
    %% Statistics
    %Label Glomerular Regions
    [labeled num]=bwlabeln(Holes, 6);
    
    labeled=single(labeled);
    if counter == 35
        labeledgloms = labeled;
    end
    
    %calculate area histogram
    num=max(labeled(:));
    num=max(num);
    areas = hist(labeled(:),num+1);
    
    areas(1) = 0;
    area_max = 1000; %max(areas);
    area_range=(1:area_max);
    [area_hist, area_axis]=hist(areas,area_max);
    stored_hist(counter,:) = area_hist;
    %Total number of detected objects in area_range
    disp(['sum(area_hist(12:100)):  ', num2str(sum(area_hist(12:100)))])
    
    %Diameter histogram from area data assuming spherical glomeruli
    [diams x]=hist((6/pi*areas(areas>2)).^(1/3),1:0.1:15); diams(end)=0;
    stored_diams(counter,:) = diams;
    
    processingcount=processingcount + 18;
    % double(processingcount)/double(totalprocessing);
    
    if processingcount == 1140
        hProgress = progressbar( hProgress, double(processingcount)/double(totalprocessing), 'Phase 2: Completed. Saving Data.');
    else
        message = ['Phase 2: Performing Threshold Statistics (' num2str(counter) '/50)'];
        ShowStatus(['XglomRunPipeline: ' message ],handles)
        hProgress = progressbar( hProgress,double(processingcount)/double(totalprocessing) , message);
    end
end
telapsed = toc(tstart);
save(fullfile(pwd,'Cache' ,handles.DatasetID, 'output_01.mat'),'stored_hist', 'stored_diams', 'labeledgloms')

ShowStatus(['Xglom Pipeline complete for dataset ' handles.DatasetID],handles);
progressbar( hProgress, -1 );
if ~hProgress
    disp('Progressbar removed.')
end

info = GetDICOMInfo(DatasetID);

PlotXGlomOutputDist(stored_hist,stored_diams,handles.DatasetID,info.PixelDimensions(1),ThresholdRange);









% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close(handles.figure1)

% --- Executes on button press in checkbox1.
function checkbox1_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox1


% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close(handles.figure1)

% --- Executes when selected cell(s) is changed in uitable1.
function uitable1_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to uitable1 (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)
% display eventdata
handles.DatasetSelection = eventdata.Indices;
guidata(hObject,handles);



% --- Executes during object creation, after setting all properties.
function headlinetext_CreateFcn(hObject, eventdata, handles)
% hObject    handle to headlinetext (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
handles.Xglom_version= getenv('XGLOM_VERSION'); % if not in environment variables there
                                                % was an error in Xglom.m
set(hObject,'String',sprintf('Xglom %s Pipeline Process',char(handles.Xglom_version)));
guidata(hObject, handles);


%% GetImage
% Searches local directory 'Cache' for dataset, otherwise gets it from
% DaRIS.  Identical to GetImage in Xglom.m but without DICOM
function imgvol = GetImage(uid,tag)

% Check for dataset in Cache
imgvol=[];
[status lsoutput] =system(['ls ./Cache/' uid]);
%     if status ~= 0
%         disp(['Error no Cache directory for ' uid])
%         return
%     end
currentFolder = pwd;
addpath(fullfile(currentFolder, '/plugins/dicom_toolbox_version_2e'));

if ~isempty(strfind(tag,'depth'))
    disp(['Getting processed depth ' uid ' from local directory Cache .']);
    [status lsoutput] =system(['ls ./Cache/' uid '/depth*.mat']);
    lsoutput = regexprep(lsoutput,'\n','');
    lsoutput = regexprep(lsoutput,'\n','');
    if status == 0
        display 'Processed depth image in Cache. Loading in.'
        load(lsoutput); imgvol= depth;
        size(imgvol)
        imgvol = single(imgvol);
        imgvol = squeeze(imgvol);
    else
        display 'Error in loading processed images.'
        return %1.;
    end
    
elseif  ~isempty(strfind(tag,'output'))
    disp(['Getting output stats of ' uid ' from local directory Cache.']);
    [status lsoutput] = system(['ls ./Cache/' uid '/output*.mat']);
    lsoutput = regexprep(lsoutput,'\n','');
    lsoutput = regexprep(lsoutput,'\n','');
    if status == 0
        display 'Processed output in Cache. Loading in.'
        load(lsoutput); imgvol=labeledgloms>0;
    else
        display 'Error in loading output stats.'
        return %1.;
    end
    
end

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

ShowStatus('XglomRunPipeline: Save as CVS  - still in testing.',handles);
display(' Save as CVS  - still in testing');

img = GetImage(handles.DatasetID,'depth');
img= NormaliseImage(img);
if numel(img)<2
    ShowStatus('XglomRunPipeline: Error reading depth image.',handles);
end

%        load(['./Cache/' handles.DatasetID '/output_01.mat'])
uiwait(msgbox({'You must create a base filename to output the processed data. The standard output will be called ''<basename>-standard.csv'' and will contain columns ID, Sizes, Estimated Diamter, X-position, Y-position, and Z-position. If you create a weighted output it will be called ''<basename>-weighted.csv'' and will contain additional column data for maximum and mean pixel intensity, and the XYZ centroids will use the the raw data rather than a binary image to calculate the centre of the label glomeruli. ','',' Select OK to create a base filename.'},'Save Stats as CSV', 'modal'));

[filename, pathname] = uiputfile('*.csv', 'Select a basename with CSV extension.');
if isequal(filename,0) || isequal(pathname,0)
    ShowStatus('User pressed cancel in Save Stats',handles)
    return
else
    ShowStatus(['User selected ' fullfile(pathname, filename)],handles)
end

%% This threshold is arbitrary, it needs to be the optimal
img = NormaliseImage(img);
Holes=img>0.28;

% %% Discard single voxels
% tic; ConArray=conndef(3,'minimal');
% ConArray(2,2,2)=0;
% Holes=imdilate(Holes,ConArray).*Holes;
% Holes=single(Holes);
% toc;

%% Statistics, Label Glomerular Regions
tic;CC = bwconncomp(Holes,6);toc;

%        labeledgloms =  labelmatrix(CC);


% calculate area histogram
%    You can calculate the following properties on N-D inputs: 'Area', 'BoundingBox', 'Centroid', 'FilledArea', 'FilledImage', 'Image', 'PixelIdxList', 'PixelList', and 'SubarrayIdx
tic;
s = regionprops(CC, 'Area', 'Centroid'); % ,'BoundingBox','SubarrayIdx');
toc;
Areas = cat(1,s.Area);

% grab the centroids, 'centroid'
%s  = regionprops(CC, 'centroid');
Centroids = cat(1, s.Centroid);

%    BoundingBox = cat(1,s.BoundingBox);
%    SubarrayIdx = cat(1,s.SubarrayIdx);
hProgress  = progressbar([],0,'Please wait, writing standard labels file.' );

ShowStatus('Saving statistics of labelled image.',handles);
dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-standard.csv')),...
    'ID,Size,EstDiameter,Xpos,YPos,ZPos','delimiter','');
hProgress  = progressbar(hProgress,0.2,'Please wait, writing standard labels file.' );
dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-standard.csv')),...
    [(1:numel(Areas))' Areas (6/pi*Areas).^(1/3) Centroids],'-append');
hProgress  = progressbar(hProgress,0.5,'Please wait, writing standard labels file.' );

progressbar( hProgress, -1 );
if ~hProgress
    disp('Progressbar removed.')
end


ShowStatus('Statistics of labelled image file written.',handles);



% Ask user if they want the weighted centroids
choice = questdlg('Would you like to save the weighted centroids (this gives more accurate co-ordinates for the centre of glomeruli)?','Yes','No');
% Handle response
switch choice
    case 'Yes'
        hProgress  = progressbar([],0,'Please wait, writing weighted labels file.' );
        
        ShowStatus('Calculating weighted centroids of labelled image.',handles);
        tic;s = regionprops(CC,img.*Holes,'WeightedCentroid', 'MaxIntensity', 'MeanIntensity');toc;
        hProgress  = progressbar(hProgress,0.2,'Please wait, writing weighted labels file.' );
        
        WCentroids = cat(1,s.WeightedCentroid);
        MaxI = cat(1,s.MaxIntensity);
        MeanI = cat(1,s.MeanIntensity);
        ShowStatus('Writing weighted centroids.',handles);
        dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-weighted.csv')),...
            'ID,Size,EstDiameter,MaxIntensity,MeanIntensity,Xpos,YPos,ZPos','delimiter','');
        hProgress  = progressbar(hProgress,0.2,'Please wait, writing weighted labels file.' );
        
        dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-weighted.csv')),...
            [(1:numel(Areas))' Areas (6/pi*Areas).^(1/3) MaxI MeanI WCentroids],'-append');
        hProgress  = progressbar(hProgress,0.4,'Please wait, writing weighted labels file.' );
        progressbar( hProgress, -1 );
        if ~hProgress
            disp('Progressbar removed.')
        end
        
        ShowStatus('Weighted centroids file written.',handles);
    case 'No'
        ShowStatus('Not writing weighted centroids.',handles)
end

% Disabling labelled image
%dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-labeled-image.csv')),Areas);
% Ask user
choice = questdlg('Would you like to save the pixel index list?','Yes','No');
% Handle response
switch choice
    case 'Yes'
        ShowStatus('Calculating pixelvalues of labelled image.',handles);
        tic;s = regionprops(CC,img.*Holes,'PixelList', 'PixelValues');toc;
        
        ShowStatus('Writing pixel list.',handles);
        hProgress  = progressbar([],0,'Please wait, writing pixel list file.' );
        dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-pixellist.csv')),...
            'PixelValue,Xpos,Ypos,Zpos','delimiter','');
        dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-pixellist.csv')),...
            [s(1).PixelValues s(1).PixelList],'-append');
        for ii=2:numel(s)
            dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-pixellist.csv')),...
                [],'-append');
            dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-pixellist.csv')),...
                [s(ii).PixelValues s(ii).PixelList],'-append');
            hProgress = progressbar( hProgress, 1.0/double(numel(s)),'Please wait, writing pixel list file.');
            
        end
        progressbar( hProgress, -1 );
        if ~hProgress
            disp('Progressbar removed.')
        end
        ShowStatus('Pixel index list file written.',handles);
        
    case 'No'
        ShowStatus('Not writing pixel index list.',handles)
end

function loadPreview(handles,img)
if ~isempty(img)
    iH = imagesc(img, 'parent', handles.pipelineaxes);
    set(iH, 'hittest', 'off');colormap(gray);
    axis(handles.pipelineaxes, 'equal','tight');
    set(handles.pipelineaxes, ...
        'box'             , 'on', ...
        'xtickmode'       , 'auto', ...
        'ytickmode'           , 'auto', ...
        'interruptible'   , 'off', ...
        'cameraviewanglemode' , 'auto', ...
        'dataaspectratiomode' , 'auto', ...
        'busyaction'      , 'queue');
end


function XglomRunPipeline_version2(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
%   - Michael Eager,   (michael.eager@monash.edu)
%   - (c) 2012, Monash Biomedical Imaging, Monash University, Australia

 
handles.debug = get(handles.checkbox2,'Value');
ShowStatus(['Running Xglom (' handles.Xglom_version ') Pipeline  ... '],handles);
tstart=tic;
ShowStatus(['Xglom (' handles.Xglom_version ') Pipeline : Loading DICOM image'],handles);
splitkidney=0;

% MR pipeline variables
Xglom_MR_pipeline.version=handles.Xglom_version;
Xglom_MR_pipeline.ProcessingDate=datestr(now);
Xglom_MR_pipeline.CorrMask_threshold=1.75;
Xglom_MR_pipeline.CorrMask_sigma=3.0/60.0;
Xglom_MR_pipeline.RegularizeMRI_factor=1;  %default
Xglom_MR_pipeline.RegularizeMRI_filter_sd=2.0; %default


if isempty(handles.Img) 
    handles.Img = GetDICOMImage(handles.DatasetID);
else
    disp('XglomRunPipeline_version2 using interal image.')
end
midaxis = floor(size(handles.Img,1)/2);

addpath(fullfile(pwd,'plugins','progressbar'));


gui_active(1);      % will add an abort button
hDummy=figure;
hProgress=figure;
hProgress_0 = hProgress;
hProgress  = progressbar([],0,' Phase 1 of 2: Image Cleaning Procedure (8 steps)' );
ShowStatus('Xglom Pipeline : Phase 1 of 2: Image Cleaning Procedure (8 steps)',handles);
incr_progress = double(1)/double(60); %estimate 240 secs for imfill, 50 reps for bwconncomp
%disable original hProgress figure
close(hProgress_0);
close(hDummy)




%% Phase 1a: Grab grayscale mask
ShowStatus('Phase 1a: Grab bias mask',handles);
% Use default args for RegularizeMRI.  To set slice range  499 to 504, with
% factor = 1 use mask=RegularizeMRI(img,499,504,1)
mask = RegularizeMRI(handles.Img);

[ymask,xmask] = hist(mask(:),100);
% estimate split kidney if more than 75% of half (0.375) of the masked
% image is at the centre point
splitkidney_upperlimit = (size(mask,1) * size(mask,2)) * 0.375;
if max(ymask) > splitkidney_upperlimit
    ShowStatus('Phase 1a: Mask shows that a split kidney is being used.',handles);
    splitkidney=1;
    bmask = handles.Img;
    zbmaskind =find(bmask~=0);
    bmask(zbmaskind)=1;
end

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
loadPreview(handles,mask);
if (handles.debug == 1 )
    disp('In debug mode')
    fig=figure;imshow(mask); drawnow;
    %pause;
    disp('calling uiwait')
    uiwait(handles.figure1);
    disp(' uiwait released ')
    if ishandle(fig)
        close(fig);
    end
end
if ~hProgress
    return
end

%% Phase 1b: Apply grayscale mask to volume
ShowStatus('Phase 1b: Apply bias mask to volume',handles);
hProgress = progressbar( hProgress, incr_progress, 'Phase 1: Step 2 of 8 - Grayscale mask');

if splitkidney == 0
    img_clean = MaskVolume(handles.Img,mask);
else
    ButtonName = questdlg('Split kidney. Do you want to use the mask to correct the image? Note, the masking method prefers the whole image to be full.', ...
        'Split kidney masking question', ...
        'Yes', 'No','Cancel','No'); % No is the default
    switch ButtonName,
        case 'Yes',
            ShowStatus('Phase 1b: Appling bias mask to split kidney image',handles);
            img_clean = MaskVolume(handles.Img,mask);
        case 'Cancel',
            ShowStatus('Phase 1b: Not using the bias mask on the split kidney image',handles);
            img_clean = handles.Img;
        case 'No',
            ShowStatus('Phase 1b: Not using the bias mask on the split kidney image',handles);
            img_clean = handles.Img;
    end % switch
end

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


loadPreview(handles,squeeze(img_clean(midaxis,:,:)));
if handles.debug==1
    fig=figure;
    imagesc(squeeze(img_clean(midaxis,:,:)));
    colormap(flipud(gray));
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end
if ~hProgress
    return
end
%% Phase 1c: Apply Laplacian of Gaussian filter to cleaned image
ShowStatus('Phase 1c: Apply Laplacian of Gaussian filter to cleaned image',handles);
hProgress = progressbar( hProgress,incr_progress, 'Phase 1: Step 3 of 8 - L-o-G filter');

h = fspecial3('log',5);
Ilog = imfilter(img_clean,h,'replicate');

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

loadPreview(handles,squeeze(Ilog(midaxis,:,:)));
%imshow(squeeze(Ilog(midaxis,:,:)))
if handles.debug==1
    fig=figure;
    imagesc(squeeze(Ilog(midaxis,:,:)));
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end


%% Phase 1d:  Smooth out Ilog greater than zero with inhomogenityCorrection
ShowStatus('Phase 1:  Step 3 of 8 - Smooth out Laplace-of-Gaussian image.',handles);
hProgress = progressbar( hProgress,incr_progress, 'Phase 1:  Step 3 of 8 - Smooth out L-o-G');
tic;corrMask=inhomogenityCorrection((Ilog.^2).*(Ilog>0),3./60);toc

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

loadPreview(handles,squeeze(corrMask(midaxis,:,:)));
if handles.debug ==1
    fig=figure;
    imagesc(squeeze(corrMask(midaxis,:,:)));
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end


if ~hProgress
    return
end
%% Phase 1e: Set a mask in smoothed-out Ilog and find connected components
ShowStatus('Phase 1:  Step 4 of 8 - Find glomerular regions',handles);
hProgress = progressbar( hProgress,incr_progress, 'Phase 1:  Step 4 of 8 - Find glomerular regions');
CC = bwconncomp(corrMask>1.75);  % threshold is selected from histogram
numPixels = cellfun(@numel,CC.PixelIdxList);
[x idx]=sort(numPixels,'descend');
idx_lim=1;



if splitkidney == 0
    count_mask_regions=0;
    while  (sum(x(1:idx_lim)) < 15000000) && count_mask_regions <= 4
        % 15000000 = 0.111*512^3
        ShowStatus('XglomRunPipeline 2: Binary labels not big enough to encompass two kidney glomerular regions, increasing regions.',handles)
        idx_lim=idx_lim+1;
        count_mask_regions=count_mask_regions+1;
    end
else
    ShowStatus('XglomRunPipeline 2: Binary labels in split kidney using just the biggest label.',handles)
end

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


if ~hProgress
    return
end
%% Phase 1f:  Create binary mask
ShowStatus('Phase 1:  Step 5 of 8 - Create binary mask',handles);
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 5 of 8 - Create binary mask');
BW = ismember(labelmatrix(CC), idx(1:idx_lim));

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

loadPreview(handles,squeeze(BW(midaxis,:,:)));
%  9.016970 seconds.
%imshow(squeeze(BW(midaxis,:,:)))
if handles.debug==1
    fig=figure;imagesc(squeeze(BW(midaxis,:,:)));
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end

if ~hProgress
    return
end


%% Phase 1g: Image filling
ShowStatus('Phase 1:  Step 6 of 8 - Flood-filling algorithm.',handles);
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 6 of 8 - Flood-filling');
telapsed = toc(tstart);
imf=single(imfill(img_clean.*BW,'holes')); %default:conn=6
depth=(imf - img_clean).*BW;

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


loadPreview(handles,squeeze(depth(midaxis,:,:)));
%imshow(squeeze(depth(midaxis,:,:)))
if handles.debug==1
    fig=figure;imagesc(squeeze(depth(midaxis,:,:)));
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end


if ~hProgress
    return
end



%% Phase 1h: Normalise the depth image and save processed data
ShowStatus('Phase 1:  Step 7 of 8 - Normalising depth image',handles);
hProgress = progressbar( hProgress, incr_progress, 'Phase 1:  Step 7 of 8 - Normalising depth image');
depth = NormaliseImage(depth);
ShowStatus('XglomRunPipeline 2:  Writing depth image to local cache.',handles);
system(['[ -f ' fullfile(pwd,'Cache',handles.DatasetID,'depth_mask_01.mat') ' ] && chmod ug+w ' fullfile(pwd,'Cache',handles.DatasetID,'depth_mask_01.mat')]);

loadPreview(handles,squeeze(depth(midaxis,:,:)));
save(fullfile(pwd,'Cache',handles.DatasetID,'depth_mask_01.mat'),'depth','mask','Xglom_MR_pipeline');
system(['chmod ug+w ' fullfile(pwd,'Cache',handles.DatasetID,'depth_mask_01.mat') ]);

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


if ~hProgress
    return
end

%% Phase 2: Thresholding
ShowStatus('Phase 2: Performing Threshold Statistics (50 steps)',handles);
hProgress = progressbar( hProgress, incr_progress, 'Phase 2: Performing Threshold Statistics (50 steps)');

telapsed = toc(tstart);

ThresholdRange= 0.11:0.01:0.5;
diams_range=1:0.1:15;
stored_hist=zeros(0,1000);
stored_diams = zeros(0,141);
counter=1;
for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %% Thresholding
    Holes=depth>ThresholdRange(counter);
    
    %% Discard single voxels
    % ConArray=conndef(3,'minimal');
    % ConArray(2,2,2)=0;
    % Holes=imdilate(Holes,ConArray).*Holes;
    % Holes=single(Holes);
    
    %% Statistics, Label Glomerular Regions
    CC = bwconncomp(Holes,6);
    %    if counter == 20
    %        labeledgloms =  labelmatrix(CC);
    %    end
    
    % calculate area histogram
    areas = cellfun(@numel,CC.PixelIdxList);
    area_max = 1000; %max(areas);
    area_range=(1:area_max);
    [area_hist, area_axis]=hist(areas,area_range);
    stored_hist(counter,:) = area_hist;
    
    
    % Total number of detected objects in area_range
    disp(['sum(area_hist(1:100)):  ', num2str(sum(area_hist(1:100)))])
    
    % Diameter histogram from area data assuming spherical glomeruli
    diams_range=1:0.1:15;
    [diams x]=hist((6/pi*areas(areas>1)).^(1/3),diams_range);
    diams(end)=0;
    stored_diams(counter,:) = diams;
    
    telapsed = toc(tstart);
    disp(['Processing time (sec): ' num2str(telapsed)]);
    
    
    if ~hProgress
        return
    end
    if counter == length(ThresholdRange)
        ShowStatus('Phase 2: Completed Xglom pipeline. Saving Data.',handles);
        
        hProgress = progressbar( hProgress, incr_progress, 'Phase 2: Completed Xglom pipeline. Saving Data.');
    else
        ShowStatus(['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)'],handles);
        
        hProgress = progressbar( hProgress,incr_progress , ['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)']);
    end
end
telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

ShowStatus('Calculating optimal threshold for detecting gloms ...',handles);


sumdiams=sum(stored_diams,2);
diam_minmax = extrems(sumdiams);

if  ~isempty(diam_minmax.maxx)
    max_i = diam_minmax.maxx(1);
    max_d = diam_minmax.maxy(1);
else
    [max_d max_i] = max(sumdiams(5:end));
end
optimal_threshold = ThresholdRange(max_i);

ShowStatus(['Optimal threshold: ' num2str(optimal_threshold) '  Number of glomeruli: ' num2str(max_d)], handles);


%% Re-reun the connected components routine to get the labelled
%  excluding the single voxels
Holes = (depth> optimal_threshold);
ConArray=conndef(3,'minimal');
ConArray(2,2,2)=0;
Holes=imdilate(Holes,ConArray).*Holes;

CC = bwconncomp(Holes,6);
labeledgloms =  labelmatrix(CC);



%% Write output data to file - ensure writability by users and group

outfile= fullfile(pwd,'Cache', handles.DatasetID, 'output_01.mat');
if exist(outfile,'file')
    system(['chmod ug+w ' outfile]);
end
save(outfile,'stored_hist', 'stored_diams', 'labeledgloms','ThresholdRange', 'diams_range', 'area_range','optimal_threshold')
system(['chmod ug+w ' outfile]);

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

ShowStatus(['Xglom Pipeline complete for dataset ' handles.DatasetID],handles);
hProgress = progressbar( hProgress, -1 , 'Image Processing Complete. Processing output plots.');

%% Gather info and plot the results

info = GetDICOMInfo(handles.DatasetID);

if exist('info.PixelDimensions')
   PixelDimension = info.PixelDimensions(1);
elseif exist('info.PixelSpacing')
    PixelDimension = info.PixelSpacing(1);
elseif exist('info.SliceThickness')
    PixelDimension = info.SliceThickness;
elseif exist('info.SpacingBetweenSlices')
    PixelDimension = info.SpacingBetweenSlices;
else
  PixelDimension=0.1; %100 micron
end
    
PlotXGlomOutputDist(stored_hist,stored_diams,handles.DatasetID,...
    PixelDimension,ThresholdRange);
handles.stored_hist = stored_hist;
handles.stored_diams = stored_diams;



function XglomRunPipeline_CTversion(hObject, eventdata, handles)
% CT pipeline of Xglom
%   The pipeline takes a 3D image (assuming the image is in the local cache)
%   under its Daris ID.  CT images come in DICOM or TIFF stacks.
%   Cleaning the CT image involves gathering a histogram of image values.
%   Extrema in the histogram indicate transitions (minima) or concentrations
%   of materials in the 3D image. The CT TIFF stack has empty space, tube,
%   embedding material, un-infused tissue, mildly infused tissue, 
%   strongly infused tissue.
%
%
%  For testing purposes use:
%  >>  [img info] = GetCTtiffImage('./Cache/1008.2.18.7.1.1.2')
%  >> CTPipeline(img,info)
%
%   - Michael Eager,   (michael.eager@monash.edu)
%   - (c) 2012, Monash Biomedical Imaging, Monash University, Australia


handles.debug = get(handles.checkbox2,'Value');
ShowStatus(['Running Xglom (' handles.Xglom_version ') CT Pipeline version 1.0 ... '],handles);
tstart=tic;
ShowStatus(['Xglom (' handles.Xglom_version ')) CT Pipeline version 1.0: Load and normalise image'],handles);

if isempty(handles.Img)
    [handles.Img info ]= GetCTtiffImage(fullfile(pwd,'Cache',handles.DatasetID));
else
    info = imfinfo(fullfile(pwd,'Cache',handles.DatasetID,'0001_01_01','0001.tif'));
end

% Get DaRIS information stored in info.txt
data_cell_string = fileread(fullfile(pwd,'Cache',handles.DatasetID,'info.txt'));
data_cell_string = regexp(data_cell_string,'\s+','split');

for i=1:length(data_cell_string)
    if ~isempty(strfind(data_cell_string{i},'NZ='))
        handles.NumSlides =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'DX='))
        handles.PixelResolution =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'DZ='))
        handles.SliceSeperation =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'TX='))
        handles.XTiles =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'TY='))
        handles.YTiles =  str2double(data_cell_string{i}(4:end));
    end
end

nimg=NormaliseImage(single(handles.Img));

midZaxis = floor(size(nimg,3)/2);



gui_active(1);      % will add an abort button
hDummy=figure;
hProgress=figure;
hProgress_0 = hProgress;
hProgress  = progressbar([],0,'CT Pipeline Phase 1 of 2: Image Cleaning Procedure (8 steps)' );
ShowStatus('Xglom CT Pipeline version 1.0: Phase 1 of 2: Image Cleaning Procedure (8 steps)',handles);
incr_progress = 0.01; %double(1)/double(60); %estimate 240 secs for imfill, 50 reps for bwconncomp
%disable original hProgress figure
close(hProgress_0);
close(hDummy)


%% Phase 1a: Grab grayscale mask
ShowStatus('CT Pipeline Phase 1 Step 1 of 6: Measure image statistics',handles);


% % calculate a histogram of the image values
% midslice=squeeze(nimg(:,:,midZaxis-50:midZaxis+50));
% %reduce_img = reducevolume(nimg,[10 10 10]);
% [hOrig xorig] = hist(midslice(:),300);
% clear midslice;
% % smooth the histogram
% hsmooth = smooth(hOrig,20);

% Faster histogram method with reduced volume
Vreduced = reducevolume(nimg,[10 10 10]);
tic;[hOrig xorig] = hist(Vreduced(:),300);toc
hsmooth = smooth(hOrig,20);


% find the maxima and minima
out = extrems(hsmooth);

minextrems = out.minx( find( (xorig(out.minx) > 0.05) &   (xorig(out.minx) < 0.5) ));
maxextrems = out.maxx( find( (xorig(out.maxx) > 0.1) &   (xorig(out.maxx) < 0.5) ));

if numel(minextrems) >= 2 && numel(maxextrems) >= 2
    
    
    % find the two minima indicating the tube boundary
    % tube_boundary = x(out.minx( find( (x(out.minx) > 0.25) &   (x(out.minx) < 0.5) )))
    tube_boundary = xorig(minextrems(1:2));
    tube = xorig(maxextrems(1));
    tissue = xorig(maxextrems(2));
    
else
    error('CT pipeline unable to calculate reasonable image ranges in histogram');
end


telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


if (handles.debug == 1 )
    % x, log10(h),'-b',
    hHist=figure;
    plot(xorig,log10(hsmooth),'-k', 'DisplayName','Smoothed');
    hold('on')
    plot(xorig,log10(hOrig),'--k','linewidth',0.5,'DisplayName','Original');
    plot(xorig(out.minx),log10(out.miny),'*r','Markersize',10,'DisplayName',[])
    plot(xorig(out.maxx),log10(out.maxy),'*g','Markersize',10,'DisplayName',[])
    hold('off'); legend('show')
    title('Histogram of raw image (with min/max points)')
    fig=figure;
    imshow(nimg(:,:,midZaxis));colormap(gray);axis image
    title('Normalised raw image')
    %pause;
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
    if ishandle(hHist)
        close(hHist);
    end
    
    %   ShowImage(maske,handles);
    %   Continue();
end
if ~hProgress
    return
end

clear hOrig

%% Phase 1b:
ShowStatus('CT Pipeline Phase 1 Step 2 of 6: Finding embedding solution region. (Est. time 6 minutes.)',handles);
hProgress = progressbar( hProgress, incr_progress*4, 'CT Pipeline Phase 1 Step 2 of 6: Finding embedding solution region. (Est. time 6 minutes.)');


CC = bwconncomp( (nimg>tube_boundary(2)));
numPixels = cellfun(@numel,CC.PixelIdxList);
[xPix idx]=sort(numPixels,'descend');

%Elapsed time is 298.715674 seconds.


telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

if handles.debug
    fig=figure;
    imagesc((nimg(:,:,midZaxis)>tube_boundary(2))); colormap(gray);axis image
    title('Thresholded Image for Connected Component Analysis')
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end


if ~hProgress
    return
end

%% Phase 1c: Fit bounding box around region and create cropped image
ShowStatus('CT Pipeline Phase 1 Step 3 of 6: Bounding box of the biggest region',handles);
hProgress = progressbar( hProgress,incr_progress*40, 'CT Pipeline Phase 1 Step 3 of 6: Get the bounding box of the biggest region (i.e. contents of tube). (Est. time 40 sec) ');


% Get the bounding box of the biggest region -> contents of tube
s = regionprops(CC, 'BoundingBox');
BB = s(idx(1)).BoundingBox;

% Elapsed time is 37.830972 seconds.
clear CC s


% crop the image and normalise values
cropped_img = NormaliseImage(single(nimg(floor(BB(1)):ceil(BB(1)+BB(4)),...
    floor(BB(2)):ceil(BB(2)+BB(5)),floor(BB(3)):ceil(BB(3)+BB(6)))));
clear nimg
midZaxis = floor(size(cropped_img,3)/2);

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

%imshow(squeeze(Ilog(midaxis,:,:)))
if handles.debug
    fig=figure;
    imagesc(squeeze(cropped_img(:,:,midZaxis)));colormap(gray);axis image
    title('Cropped image: Bounding Box from biggest connected component')
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end


%% Phase 1d:  Smooth out Ilog greater than zero with inhomogenityCorrection
ShowStatus('CT Pipeline Phase 1 Step 4 of 6: Crop the CT image.',handles);

hProgress = progressbar( hProgress,incr_progress*3, 'CT Pipeline Phase 1 Step 4 of 6: Crop the CT image.');


crop_range_img = cropped_img;
crop_range_img(cropped_img < tissue) = tissue;
% crop_range_img(cropped_img>0.6) = 0.6;
crop_range_nimg = NormaliseImage(crop_range_img);

clear crop_range_img
midZaxis = floor(size(crop_range_nimg,3)/2);


%tic;corrMask=inhomogenityCorrection((Ilog.^2).*(Ilog>0),3./60);toc

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

%imshow(squeeze(corrMask(midaxis,:,:)))
if handles.debug ==1
    fig=figure;
    imagesc(crop_range_nimg(:,:,midZaxis));colormap(gray);axis image
    title('Scaled and normalised cropped image')
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end
if ~hProgress
    return
end



%% Phase 1g: Image filling
ShowStatus('CT Pipeline Phase 1 Step 5 of 6: Laplace-of-Gassian filter.',handles);
hProgress = progressbar( hProgress, incr_progress, 'CT Pipeline Phase 1 Step 5 of 6: Laplace-of-Gassian filter. (Est. time 5 minutes.)');

% create image filter of laplacian-of-gaussian operator
hLoG = fspecial3('log',5);
Ilog_tissue= imfilter((crop_range_nimg).*(cropped_img>(tissue)),hLoG,'replicate');

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
%
% %imshow(squeeze(depth(midaxis,:,:)))
if handles.debug==1
    fig=figure;
    imagesc(squeeze(Ilog_tissue(:,:,midZaxis)));%colormap(gray);
    axis image
    title('Laplace-of-Gaussian (Blob detector filter) image')
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end
if ~hProgress
    return
end



%% Phase 1h: Normalise the depth image and save processed data
ShowStatus('CT Pipeline  Phase 1  Step 6 of 6: Normalise LoG output and Save the processed image',handles);
hProgress = progressbar( hProgress, incr_progress*19, 'CT Pipeline  Phase 1  Step 6 of 6: Save the processed image');

Ilog_tmp = (Ilog_tissue.*(Ilog_tissue<0).*(cropped_img>(tissue*1.01)));
Ilog_tmp(Ilog_tissue < -0.05)=-0.05;

depth = 1- NormaliseImage(Ilog_tmp);

%depth = crop_range_nimg;  % NormaliseImage(depth);

if handles.debug==1
    fig=figure;
    
    imagesc(squeeze(depth(:,:,midZaxis)));colormap(gray);axis image ;
    title('Depth image (Final processed image before threshold phase)');
    uiwait(handles.figure1);
    if ishandle(fig)
        close(fig);
    end
end
if ~hProgress
    return
end



telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
if exist(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'file')
    system(['chmod ug+w ' fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat')]);
end
save(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'depth','crop_range_nimg','tissue', 'tube_boundary','-v7.3'); %, 'Ilog_tissue'
if exist(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'file')
    system(['chmod ug+w ' fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat')]);
end
clear crop_range_nimg cropped_img Ilog_tissue

if ~hProgress
    return
end

%% Phase 2: Thresholding
ShowStatus('Phase 2: Performing Threshold Statistics (50 steps)',handles);
hProgress = progressbar( hProgress, incr_progress/10, 'Phase 2: Performing Threshold Statistics (50 steps)');

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
ThresholdRange= 0.11:0.01:0.4;
stored_hist=zeros(0,1000);
diams_range=1:0.1:300;
stored_diams = zeros(0,length(diams_range));
area_max = 1000; %max(areas);
area_range=(1:area_max);

for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %% Thresholding
    Holes=depth>ThresholdRange(counter);
    
    %% Discard single voxels
    % ConArray=conndef(3,'minimal');
    % ConArray(2,2,2)=0;
    % Holes=imdilate(Holes,ConArray).*Holes;
    % Holes=single(Holes);
    
    %% Statistics, Label Glomerular Regions
    CC = bwconncomp(Holes,6);
    %    if counter == 20
    %        labeledgloms =  labelmatrix(CC);
    %    end
    
    % calculate area histogram
    areas = cellfun(@numel,CC.PixelIdxList);
    
    area_hist=hist(areas,area_range);
    stored_hist(counter,:) = area_hist;
    
    
    % Total number of detected objects in area_range
    disp(['sum(area_hist(1:100)):  ', num2str(sum(area_hist(5:200)))]);
    
    % Diameter histogram from area data assuming spherical glomeruli
    diams=hist((6/pi*areas).^(1/3),diams_range); diams(end)=0;
    stored_diams(counter,:) = diams;
    
    %disp(['sum(diams(1:100)):  ', num2str(sum(diams(1:100)))]);
    
    
    if ~hProgress
        return
    end
    if counter == length(ThresholdRange)
        ShowStatus('Phase 2: Completed Xglom pipeline. Saving Data.',handles);
        
        hProgress = progressbar( hProgress, incr_progress, 'Phase 2: Completed Xglom pipeline. Saving Data.');
    else
        ShowStatus(['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)'],handles);
        
        hProgress = progressbar( hProgress,incr_progress , ['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)']);
    end
end
telapsed = toc(tstart);

ShowStatus('Calculating optimal threshold for detecting gloms ...',handles);

% % Find the best threshold assuming the gloms are around 100 voxels in size
% areas_estimate = sum(stored_hist(:,100:150)');
% Areas_minmax= extrems(areas_estimate);
%
% % h100= smooth(areas100,3);
% % diams_tmp = stored_diams';
% %
% % sumdiams=sum(stored_diams');
% % diam_minmax = extrems(sumdiams);
% %
% % areas_tmp = stored_hist';
% % areas_tmp(end,:) = 0;
% % areas_tmp(1:20,:) = 0;
% % hSmooth = smooth(areas_tmp(5:300,11),20);
% % hSmooth = smooth(hSmooth,20);
% % sumhist=sum(areas_tmp);
% % Areas_minmax = extrems(hSmooth)
% if  ~isempty(Areas_minmax.maxx)
%     max_i = Areas_minmax.maxx(1);
%     max_d = Areas_minmax.maxy(1);
% else
%     [max_d max_i] = max(areas_estimate);
% end
% optimal_threshold = ThresholdRange(max_i);
%
% % Find the mean and distribution
% % this process should be done using a fitting method
% %   opt_hist = stored_hist(max_i,:)
% %   hSmooth = smooth(opt_hist,10,'lowess',10);hSmooth = smooth(hSmooth,10);
% %   GlomAreas_minmax = extrems(hSmooth)
% %
% % if  ~isempty(GlomAreas_minmax.maxx) && ~isempty(GlomAreas_minmax.minx)
% %
% %
% %    GlomOptimum_mean = GlomAreas_minmax.maxx(1);
% %    GlomOptimim_sd = GlomAreas_minmax.maxy(1);
% %
% %    GlomOptimum_sum = stored_hist(max_i,GlomAreas_minmax.minx(1): 300)
% %
% % else
% %     [max_d max_i] = max(areas_estimate);
% % end
%
% opt_hist = stored_hist(max_i,:);
% hSmooth = smooth(opt_hist,10,'lowess',10);
% x=1:numel(hSmooth);
% fvoxels = fit(x(:),hSmooth(:),'gauss3')
% if fvoxels.b2 < 300 && fvoxels.b2 > 40
%     GlomVOptimum.mean = fvoxels.b2;
%     GlomVOptimum.sd = sqrt((fvoxels.c2)/2);
% elseif  fvoxels.b1 < 300 && fvoxels.b1 > 40
%     GlomVOptimum.mean = fvoxels.b1;
%     GlomVOptimum.sd = sqrt((fvoxels.c1)/2);
% elseif  fvoxels.b3 < 300 && fvoxels.b3 > 40
%     GlomVOptimum.mean = fvoxels.b3;
%     GlomVOptimum.sd = sqrt((fvoxels.c3)/2);
% else
%     disp('Fitting procedure failed to find glom size range between 300 and 40 voxels')
%     GlomVOptimum.mean = 0;
%     GlomVOptimum.sd = 0;
% end
% GlomVOptimum.total = sum(hSmooth(24:250))
% disp(GlomVOptimum)
% voxel_fit_output = evalc('disp(fvoxels)');
%
% x=diams_range(300:1500);%
% opt_diams = stored_diams(max_i,300:1500)
% hSmooth = smooth(opt_diams,10,'lowess');
% fdiams = fit(x(:),hSmooth(:),'gauss3')
% if fdiams.b2 < 200 && fdiams.b2 > 40
%     GlomDiamOptimum.mean = fdiams.b2;
%     GlomDiamOptimum.sd = sqrt((fdiams.c2)/2);
% elseif  fdiams.b1 < 200 && fdiams.b1 > 40
%     GlomDiamOptimum.mean = fdiams.b1;
%     GlomDiamOptimum.sd = sqrt((fdiams.c1)/2);
% elseif  fdiams.b3 < 200 && fdiams.b3 > 40
%     GlomDiamOptimum.mean = fdiams.b3;
%     GlomDiamOptimum.sd = sqrt((fdiams.c3)/2);
% else
%     disp('Fitting procedure failed to find glom diam size range between 40 and 200 microns')
%     GlomDiamOptimum.mean = 0;
%     GlomDiamOptimum.sd = 0;
% end
% GlomDiamOptimum.total = sum(hSmooth(200:1000))
%
% disp(GlomDiamOptimum)
%
% diams_fit_output = evalc('disp(fdiam)');
% % if handles.debug==1
% %  plot(stored_hist(max_i,:)); xlim ([1 400])
% %   hold('on'); %plot(hSmooth,'-r')
% %   plot(GlomAreas_minmax.maxx,GlomAreas_minmax.maxy,'*r','Markersize',10)
% %   plot(GlomAreas_minmax.minx,GlomAreas_minmax.miny,'*g','Markersize',10)
% %   hold('off')
% % end

GlomStats = AnalyseThresholdStatistics(stored_diams, stored_hist, handles.PixelResolution, ThresholdRange, diams_range);


ShowStatus(['Optimal threshold: ' num2str(GlomStats.optimal_threshold) '  Number of glomeruli (Voxel Est.)  ' num2str(GlomStats.GlomVoxelOptimum.total) '(Diam Est.) ' num2str(GlomStats.GlomDiamOptimum.total)], handles);
hProgress = progressbar( hProgress, incr_progress, 'Xglom CT pipeline. Getting optimum labeled image. Saving Data.');


% re-reun the connected components routine to get the labelled
%CC = bwconncomp((depth > optimal_threshold),6);
%areas = cellfun(@numel,CC.PixelIdxList);
%idx = find ( (areas < 350) & (areas > 20));
%labelebgloms = labelmatrix(CC.PixelIdxList{idx})

BW = depth> GlomStats.optimal_threshold;
BW2 = bwareaopen(BW, 20);
BW3 = bwareaopen(BW2, 350);
BW4 = BW2 - BW3;
labeledgloms =  bwlabeln(BW4);



%% Write output data to file - ensure writability by users and group
system(['[ -f ./Cache/' handles.DatasetID '/output_01.mat ] && chmod ug+w ./Cache/' handles.DatasetID '/output_01.mat']);
save(fullfile(pwd,'Cache', handles.DatasetID, 'output_01.mat'),'stored_hist', 'stored_diams', 'labeledgloms','ThresholdRange','GlomStats','diams_range','-v7.3')
system(['chmod ug+w ./Cache/' handles.DatasetID '/output_01.mat']);



ShowStatus(['Xglom Pipeline complete for dataset ' handles.DatasetID],handles);
hProgress = progressbar( hProgress, -1 , 'Image Processing Complete. Processing output plots.');



PlotXGlomOutputDist(stored_hist,stored_diams,handles.DatasetID,handles.PixelResolution,ThresholdRange);
handles.stored_hist = stored_hist;
handles.stored_diams = stored_diams;

telapsed = toc(tstart);
disp(['Xglom CT Pipeline Processing time (sec): ' num2str(telapsed)])



% --- Executes on button press in checkbox2.
function checkbox2_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox2
handles.debug = get(hObject,'Value');

if handles.debug == 1
    set(handles.debug_continue_button,'visible','on','enable','on');
else
    set(handles.debug_continue_button,'visible','off','enable','off');
end

% --- Executes on button press in debug_continue_button.
function debug_continue_button_Callback(hObject, eventdata, handles)
% hObject    handle to debug_continue_button (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.debug = 0;
uiresume(gcbf)
guidata(hObject,handles)



function ShowStatus(message,handles)
disp(message)
set(handles.statusText,'string',message);
drawnow;
% future work will write message to log file


% --- Executes on button press in pushbutton6.
%  Save Stats as CSV
function pushbutton6_Callback(~, ~, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
ShowStatus('Loading depth image.',handles);
img = GetImage(handles.DatasetID,'output');
ShowStatus('Normalising depth image.',handles);
img = NormaliseImage(img);
flag=0;
if numel(img)<2
    flag=1;
end
ShowStatus('Loading output stats.',handles);
load(fullfile('Cache', handles.DatasetID, 'output_01.mat'));

uiwait(msgbox({'You must create a base filename to output the glomeruli size histogram and diameter histogram. The two outputs will be called ''<basename>-sizes.csv'' or ''<basename>-diams.csv''. ','',' Select OK then create a base filename.'},'Save Stats as CSV', 'modal'));
[filename, pathname] = uiputfile('*.csv', 'Select a basename with CSV extension.');
if isequal(filename,0) || isequal(pathname,0)
    ShowStatus('User pressed cancel in Save Stats',handles)
    return
else
    ShowStatus(['User selected ', fullfile(pathname, filename)],handles)
end
disp('Saving Area histograms')
dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-sizes.csv')),stored_hist)
disp('Saving Diameter histograms')
dlmwrite(fullfile(pathname,regexprep(filename,'.csv','-diams.csv')),stored_diams)


% --- Executes during object creation, after setting all properties.
function debugwindow_CreateFcn(hObject, ~, handles)
% hObject    handle to debugwindow (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate debugwindow
kidney = imread('kidney.png');
image(kidney);
%ShowImage(kidney,handles)



% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, ~, handles)
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called


kidney=imread('kidney.png');
imshow(kidney); axis image; axis off;
handles.debugwindow=hObject;
guidata(hObject,handles)


function ShowImage(img,handles)

imshow(img,'Parent',handles.debugwindow);
axis image;axis off;
drawnow;


% --- Executes during object creation, after setting all properties.
function pipelineaxes_CreateFcn(hObject, eventdata, handles) %#ok<*INUSL>
% hObject    handle to pipelineaxes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate pipelineaxes
disp('creating pipelineaxes')
kidney=imread('kidney.png');
% %handles.splashimage=kidney;
% ax = get(handles.axes1)
% set(ax,'CData',squeeze(img(256,:,:)),'axis','image');
% %  if size(img,3) ==512
% %      h=imagesc(squeeze(img(256,:,:))); axis image;
% %
% %  end
imshow(kidney); axis image; axis off;
handles.debugwindow=hObject;
guidata(hObject,handles)
